Research Article | Open Access
Yuto Endo, Jun Tanida, Makoto Naruse, Ryoichi Horisaki, "Extrapolated Speckle-Correlation Imaging", Intelligent Computing, vol. 2022, Article ID 9787098, 8 pages, 2022. https://doi.org/10.34133/2022/9787098
Extrapolated Speckle-Correlation Imaging
Imaging through scattering media is a longstanding issue in a wide range of applications, including biomedicine, security, and astronomy. Speckle-correlation imaging is promising for noninvasively seeing through scattering media by assuming shift invariance of the scattering process called the memory effect. However, the memory effect is known to be severely limited when the medium is thick. Under such a scattering condition, speckle-correlation imaging is not practical because the correlation of the speckle decays, reducing the field of view. To address this problem, we present a method for expanding the field of view of single-shot speckle-correlation imaging by extrapolating the correlation with a limited memory effect. We derive the imaging model under this scattering condition and its inversion for reconstructing the object. Our method simultaneously estimates both the object and the decay of the speckle correlation based on the gradient descent method. We numerically and experimentally demonstrate the proposed method by reconstructing point sources behind scattering media with a limited memory effect. In the demonstrations, our speckle-correlation imaging method with a minimal lensless optical setup realized a larger field of view compared with the conventional one. This study will make techniques for imaging through scattering media more practical in various fields.
Seeing through scattering media is an important research topic in optics and photonics because of its wide range of applications. For example, microscope imaging through biological tissues and telescopic observation through atmospheric turbulence are longstanding issues in biomedicine and astronomy, respectively [1–4]. Recent advancements in optics and information science have driven studies for imaging through strongly scattering media where conventional methods assuming ballistic photons are difficult to apply [5–7].
Wavefront shaping based on feedback and object retrieval with a transmission matrix are established approaches for imaging through scattering media. In the wavefront shaping approach, a scattering pattern on an image sensor inside the scattering medium is fed back to a spatial light modulator outside the medium to focus light on the sensor [8–10]. In the transmission matrix approach, the scattering process is described as a matrix, and the object is recovered from a single captured image with the inversion of the matrix [11–13]. The issues with those approaches are invasiveness and complex optical setups for the feedback to optimize the wavefront or the observation of the transmission matrix.
Speckle-correlation imaging is a promising noninvasive approach with a minimal optical setup to address the above issues [14, 15]. This approach assumes shift invariance of the scattering impulse response called the memory effect, and the object is recovered from the autocorrelation of the captured speckle pattern by using a phase retrieval algorithm [16–19]. It is also extendable to multidimensional imaging while maintaining its optical simplicity [20–22].
An issue with speckle-correlation imaging is the small field of view due to a limited range of the memory effect when the scattering medium is thick. To overcome this difficulty, ptychographic methods have been introduced, although these require multishot measurements and additional hardware components [23, 24]. Other interesting methods for single-shot imaging include decomposition of multiplexed speckle correlations and localization of the speckle correlation [25, 26]. These methods require some specific optical conditions, such as isolated objects or a near-field setting.
Here, we propose and demonstrate a method for extending the field of view of single-shot speckle-correlation imaging to address the above issues. The proposed method takes account of the decay of the speckle correlation under a limited memory effect and extrapolates the correlation in the reconstruction process. Our method is readily applicable to conventional speckle-correlation imaging methods without any optical modifications. Therefore, this study will contribute to various imaging applications where scattering causes limitations.
In our method, an object is observed as a captured image on an image sensor through a scattering process without any additional device, as shown in Figure 1. This measurement process is written with a shift-variant point spread function (PSF) as where and are two-dimensional coordinates on the sensor and object planes, respectively. Here, is the response at on the sensor plane from the impulse at on the object plane.
When the memory effect is limited, in Equation (2), which is the correlation between the two scattering PSFs from the two impulse positions and , respectively, is described as where is background noise, is a decay function of the correlation, and is a parameter for the decay function [17, 18, 27]. By using the paraxial approximation, this decay function is written as
Here, is an unknown in this study, where is the wave number, is the thickness of the scattering medium, and is the distance between the scattering medium and the sensor plane. A large represents a limited memory effect caused by a thick scattering medium.
By substituting Equation (3) into Equation (2) and introducing a variable , the autocorrelation of the captured image is rewritten as where is background noise and denotes an operator of the autocorrelation. Therefore, by ignoring the background noise, the autocorrelation of the captured image is the product of the autocorrelation of the object and the decay function .
In this study, we simultaneously calculate the estimations of both the object and the decay parameter by solving the following optimization problem: where is an error function with the forward model in Equation (5). The inverse problem of Equation (5) is ill-conditioned and is difficult to solve by itself, so we introduce a penalty function for the object to regulate the inversion. Here, we define the error function as follows: where is the support of the object and is the support of the autocorrelation. Both of the supports are binary patterns. The object’s support constrains the estimated object area . The autocorrelation’s support is introduced to remove the measurement noise on the captured image because the noise concentrates at the central peak of the autocorrelation and is removed by a central hole on . Here, the background noise can be ignored in the forward model in Equation (5) by computational removal after taking the autocorrelation of the speckle image [14, 15, 28].
The penalty function is composed of the following three subfunctions: where , , and are penalties of the nonnegativity, the upper limitation with the threshold , and the sparsity with the norm, respectively. Here, , , and are tuning parameters for each of the subfunctions in the optimization process. These subfunctions for the object’s penalties are defined as where “min” is an operator outputting the smaller values of the two parenthesized variables and “max” is an operator outputting the larger values of the two parenthesized variables, respectively. Our algorithm in this study supposes sparse objects, since we introduced the sparsity regularization in Equation (12), and is therefore not applicable to dense objects. However, the assumption of the object’s sparsity is acceptable in specific applications, such as astronomical observation and superresolution fluorescence microscopy, where objects are sparse [3, 29–31]. The limitation given by the sparsity regularization is demonstrated in the numerical analysis, and ways of mitigating it are mentioned in the conclusion.
We solve the optimization problem in Equation (6) based on an iterative gradient descent algorithm in which the decayed autocorrelation is extrapolated. The gradient descent steps for each and are written as where is the index of the iterations and and are parameters for tuning the steps.
We derive the partial derivatives in Equations (13) and (14) by using the chain rule . In Equation (13), the partial derivative in the second term on the right side is calculated as where and denote the Fourier transform and the inverse Fourier transform, respectively. We calculate the autocorrelation in the Fourier domain for the first partial derivative in Equation (15). The partial derivative in the third term is written as
Each result of the derivatives in Equations (15) and (16) does not have imaginary parts when is composed of real numbers, and thus, the realness of the updated object is satisfied. In Equation (14), the partial derivative in the second term on the right side is calculated as
For the first partial derivative in Equation (17), we differentiate Equation (4) with respect to the decay parameter . The updating processes in Equations (13) and (14) iterate until the cost function in Equation (6) converges, and then, the limited autocorrelation is extrapolated.
3. Demonstration and Discussion
We conducted a simulation to quantitatively analyze the proposed method compared with the conventional one. In this simulation, the pixel count of the object was , and point sources were randomly located on a pixel central area of the object. The density of the point sources in the area was varied from 0.005 to 0.03 at intervals of 0.005. The decay parameter was also varied from 0 to 0.14 at intervals of 0.02, where corresponded to the decay-free memory effect. The reconstruction results at each density and decay parameter were evaluated with the peak signal-to-noise ratio (PSNR) [25, 33], which was defined as
As a reference, a reconstruction algorithm without considering the decay function was also employed. In this conventional approach, the estimated decay parameter and the tuning parameter were set to 0 in the updating step of Equation (14) in the proposed algorithm so that the decay function could be ignored. We refer to this simplified version as the conventional algorithm.
In both the proposed and conventional algorithms, the object’s support was a binary square of pixels, and the autocorrelation’s support was a circular hole with a diameter of 5 pixels. The object was recovered from the decayed autocorrelation with 60,000 iterations in the algorithms. The tuning parameters in the algorithms were set to , , , and . in the proposed algorithm was set to . We empirically tuned these parameters by supposing a unimodal potential map of the inverse problem. This parameter tuning process may be automated by using a grid search to find the parameters that minimize the error function in Equation (7). The algorithms started from ten combinations of random patterns for the initial and random values for the initial , and one result that achieved the minimal error function was chosen [14, 15]. The calculation time for this reconstruction process was about 20 minutes using MATLAB with an Intel Xeon Gold 6254 processor running at 3.10 GHz and equipped with 376 GB of RAM. This calculation time may be reduced by using parallel processing with a GPU and accelerating the gradient step by employing algorithms that use momentum .
The PSNRs depending on the normalized decay parameter and the density are shown in Figure 2, where each PSNR was the average value calculated from ten randomly generated objects. Here, the normalization of the decay parameter was introduced to compare the simulation results and the experimental results below, where the pixel count was different from that in the simulation. Under most conditions, the proposed algorithm in Figure 2(a) achieved higher PSNRs compared with the conventional one in Figure 2(b). These results showed the advantage of the proposed algorithm over the conventional one. As exceptions, the performance of the conventional algorithm was better when the decay parameter was 0. This was because the conventional algorithm assumes a decay-free memory effect. In both the proposed and conventional algorithms, the PSNRs were low when the object was dense because such objects conflicted with the sparsity regularization in Equation (12) and the autocorrelation’s contrast of such objects degraded. The PSNRs in the proposed algorithm decreased when the density was too low or the decay was too large, as shown in Figure 2(a). In the former case, sampling points for the decay of the autocorrelation were too sparse and insufficient to estimate the decay because the autocorrelation of such objects was sparse. The latter case was an underdeterminant case because the width of the autocorrelation was too small to estimate the object. If we define successful reconstructions as those having a PSNR of 70 dB in Figure 2, the number of conditions achieving this PSNR with the proposed method was 13, and that with the conventional method was 3. In this case, the proposed method increased the number of successful conditions to 4.3 times that of the conventional method.
In the proposed algorithm, the decay parameter was estimated simultaneously with the reconstruction of the object . The root-mean-square errors (RMSEs) between the normalized original and estimated decay parameters under each condition in Figure 2(a) are shown in Figure 3. By comparing these two figures, the PSNR of the object estimation was high when the RMSE of the decay parameter estimation was low. This result shows that the estimations of the object and the decay parameter worked together in the proposed algorithm. Therefore, the object reconstruction was aborted if the estimation of the decay parameter failed, and vice versa.
We experimentally demonstrated the proposed method with the optical setup shown in Figure 4. Collimated light from a white light-emitting diode (LED, XLamp CXA2520 manufactured by CREE) passed through a bandpass filter (BPF, HMZ0530 manufactured by ASAHI SPECTRA, central wavelength: 530 nm, bandwidth: 10 nm) to implement spatially incoherent monochromatic illumination. A diffuser (KHYP1-12 manufactured by Optical Solutions, circular diffusion angle: 5°) was illuminated by the spatially incoherent light. An object, which was a piece of aluminum foil with fifteen holes, as shown in Figure 5(c), was illuminated by the diffused light. Light passing through the object was captured by an image sensor (DMK38UX253 manufactured by The Imaging Source, pixel count: , pixel pitch: 3.45 μm) through another diffuser (KHYP1-12 manufactured by Optical Solutions, circular diffusion angle: 5°) without any imaging optics. The decay parameter is variable depending on the distance between the object and the second diffuser, and becomes larger when the distance is shorter. In the experiment, the distance between the object and the second diffuser was 40 mm, and the distance between the second diffuser and the image sensor was 27 mm, respectively.
The captured image is shown in Figure 5(a), where the point sources are not recognizable. To compensate for the shading effect, the captured image was divided by a blurred version of the captured speckle image with low-pass filtering. Then, the autocorrelation of the compensated image was calculated. The central region of pixels of the autocorrelation was clipped after subtracting the background noise in Equation (5), which was determined with the maximal pixel value outside of the clipped area. Then, the pixel count was resized to pixels to reduce the computational cost. The autocorrelation result with the support is shown in Figure 5(b). The proposed and conventional algorithms were applied to the resultant image, as mentioned in the simulation section above. In the experiment, the object’s support was a binary square of pixels, and the autocorrelation’s support was a circular hole with a diameter of 50 pixels. The tuning parameters were set to , , , and , in the both algorithms, and in the proposed algorithm. The other conditions were the same as the above simulation.
After estimating the object and the decay parameter in the proposed and conventional algorithms, a denoising process was additionally applied to the results of both algorithms by using a modified version of the conventional algorithm, where the decay parameter was not updated because . In the modified conventional algorithm, the penalty term in Equation (12) was replaced with to remove small noise defined by . The number of iterations in this denoising process was 30,000.
The reconstruction results obtained by the proposed and the conventional algorithms are shown in Figures 5(d) and 5(e), respectively. By comparing the object image in Figure 5(c) and these results, the proposed algorithm recovered point sources in a larger area than the conventional one, where the point sources on the lower part disappeared. The errors calculated in Equation (7) for the results obtained by the proposed and conventional algorithms were and , respectively. Therefore, the proposed algorithm outperformed the conventional one.
Furthermore, the autocorrelations of the object image in Figure 5(c) and the reconstructed images in Figures 5(d) and 5(e) are shown in Figures 5(f)–5(h), respectively. By comparing these autocorrelations, the autocorrelation estimated by the proposed algorithm was broader than that by the conventional one. This result verified our concept of speckle-correlation imaging with the extended field of view by extrapolating the limited autocorrelation. In addition, both the proposed and conventional algorithms interpolated the central peak of the autocorrelation removed by the autocorrelation support shown in Figure 5(b).
In the experiment, the proposed algorithm estimated the normalized decay parameter as , and the estimated decay function is shown in Figure 6. The actual correlations of the PSFs depending on the lateral position were experimentally measured by scanning a point source with intervals of 0.1 mm. The actual normalized decay parameter was calculated to be from the experimentally observed correlations by solving based on the least-squares method. The actual decay function and the correlation plot are also shown in Figure 6. In this case, the RMSE of the normalized estimated decay parameter was 0.84. The density of the object in the experiment was . The normalized decay parameter and density in the experiment corresponded to the condition where the proposed method was applicable but the conventional method was not applicable, as shown in Figure 2. In addition, the experimental RMSE was slightly higher than the simulated one at the normalized decay parameter and the density shown in Figure 3, and this RMSE was permissible for the object reconstruction. From these results of the reconstructed object and the estimated decay parameter , it was verified that the proposed method extended the field of view under the limited memory effect.
We presented and demonstrated a method for extending the field of view of single-shot speckle-correlation imaging under a limited memory effect. We modeled the speckle correlation taking into account such a situation. Then, we derived the inverse process to the object and the decay parameter of the correlation from the speckle correlation. The object and the decay parameter were simultaneously estimated based on the gradient descent algorithm, and as a result, the limited autocorrelation was extrapolated. We verified the proposed method in both a simulation and an experiment. The simulation result showed that the proposed algorithm outperformed the conventional algorithm when the memory effect was limited. In the experiment, the proposed algorithm recovered a larger field of view compared with the conventional algorithm, and its estimated decay parameter agreed with the actually measured one.
The proposed method is readily applied to conventional speckle-correlation imaging without any optical modification. Therefore, it enhances the practicality of imaging techniques for seeing through scattering media with noninvasiveness and a minimal lensless setup. In this study, we assumed spatially sparse point sources as objects for the proof of concept, and this assumption is acceptable in several applications, such as astronomical observation and superresolution fluorescence microscopy [3, 29–31]. The sparsity assumption in our method may be mitigated by using phase retrieval based on compressive sensing and deep learning [35–38]. It is also extendable to multidimensional speckle-correlation imaging [20–22]. Therefore, the proposed method will contribute to imaging applications in various fields, such as biomedicine, astronomy, and security.
Data used in this study is available from the corresponding author upon reasonable request.
Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this article.
This work was supported by Japan Society for the Promotion of Science KAKENHI Grant Numbers JP20H02657, JP20K05361, and JP20H05890.
- V. Ntziachristos, “Going deeper than microscopy: the optical imaging frontier in biology,” Nature Methods, vol. 7, no. 8, pp. 603–614, 2010.
- N. Ji, “Adaptive optical fluorescence microscopy,” Nature Methods, vol. 14, no. 4, pp. 374–380, 2017.
- R. Davies and M. Kasper, “Adaptive optics for astronomy,” Annual Review of Astronomy and Astrophysics, vol. 50, no. 1, pp. 305–351, 2012.
- A. T. Watnik and D. F. Gardner, “Wavefront sensing in deep turbulence,” Optics & Photonics News, vol. 29, no. 10, pp. 38–45, 2018.
- A. P. Mosk, A. Lagendijk, G. Lerosey, and M. Fink, “Controlling waves in space and time for imaging and focusing in complex media,” Nature Photonics, vol. 6, no. 5, pp. 283–292, 2012.
- R. Horstmeyer, H. Ruan, and C. Yang, “Guidestar-assisted wavefront-shaping methods for focusing light into biological tissue,” Nature Photonics, vol. 9, no. 9, pp. 563–571, 2015.
- S. Yoon, M. Kim, M. Jang et al., “Deep optical imaging within complex scattering media,” Nature Reviews Physics, vol. 2, no. 3, pp. 141–158, 2020.
- I. M. Vellekoop and A. P. Mosk, “Focusing coherent light through opaque strongly scattering media,” Optics Letters, vol. 32, no. 16, pp. 2309–2311, 2007.
- I. M. Vellekoop, A. Lagendijk, and A. P. Mosk, “Exploiting disorder for perfect focusing,” Nature Photonics, vol. 4, no. 5, pp. 320–322, 2010.
- O. Katz, E. Small, Y. Bromberg, and Y. Silberberg, “Focusing and compression of ultrashort pulses through scattering media,” Nature Photonics, vol. 5, no. 6, pp. 372–377, 2011.
- S. Popoff, G. Lerosey, M. Fink, A. C. Boccara, and S. Gigan, “Image transmission through an opaque material,” Nature Communications, vol. 1, no. 1, 2010.
- A. Liutkus, D. Martina, S. Popoff et al., “Imaging with nature: compressive imaging using a multiply scattering medium,” Scientific Reports, vol. 4, no. 1, p. 5552, 2014.
- R. Horisaki, R. Takagi, and J. Tanida, “Learning-based imaging through scattering media,” Optics Express, vol. 24, no. 13, pp. 738–743, 2016.
- J. Bertolotti, E. G. V. Putten, C. Blum, A. Lagendijk, W. L. Vos, and A. P. Mosk, “Non invasive imaging through opaque scattering layers,” Nature, vol. 491, no. 7423, pp. 232–234, 2012.
- O. Katz, P. Heidmann, M. Fink, and S. Gigan, “Non-invasive single-shot imaging through scattering layers and around corners via speckle correlations,” Nature Photonics, vol. 8, no. 10, pp. 784–790, 2014.
- J. R. Fienup, “Phase retrieval algorithms: a comparison,” Applied Optics, vol. 21, no. 15, pp. 2758–2769, 1982.
- S. Feng, C. Kane, P. A. Lee, and A. D. Stone, “Correlations and fluctuations of coherent wave transmission through disordered media,” Physical Review Letters, vol. 61, no. 7, pp. 834–837, 1988.
- I. Freund, M. Rosenbluh, and S. Feng, “Memory effects in propagation of optical waves through disordered media,” Physical Review Letters, vol. 61, no. 20, pp. 2328–2331, 1988.
- J. R. Fienup, “Phase retrieval algorithms: a personal tour,” Applied Optics, vol. 52, no. 1, pp. 45–56, 2013.
- Y. Okamoto, R. Horisaki, and J. Tanida, “Noninvasive three-dimensional imaging through scattering media by three-dimensional speckle correlation,” Optics Letters, vol. 44, no. 10, pp. 2526–2529, 2019.
- R. Horisaki, Y. Okamoto, and J. Tanida, “Single-shot noninvasive three-dimensional imaging through scattering media,” Optics Letters, vol. 44, no. 16, pp. 4032–4035, 2019.
- K. Ehira, R. Horisaki, Y. Nishizaki, M. Naruse, and J. Tanida, “Spectral speckle-correlation imaging,” Applied Optics, vol. 60, no. 8, pp. 2388–2392, 2021.
- G. Li, W. Yang, H. Wang, and G. Situ, “Image transmission through scattering media using ptychographic iterative engine,” Applied Sciences, vol. 9, no. 5, p. 849, 2019.
- M. Rosenfeld, G. Weinberg, D. Doktofsky, Y. Li, L. Tian, and O. Katz, “Acousto-optic ptychography,” Optica, vol. 8, no. 6, pp. 936–943, 2021.
- X. Wang, X. Jin, J. Li, X. Lian, X. Ji, and Q. Dai, “Prior-information-free single-shot scattering imaging beyond the memory effect,” Optics Letters, vol. 44, no. 6, pp. 1423–1426, 2019.
- M. Alterman, C. Bar, I. Gkioulekas, and A. Levin, “Imaging with local speckle intensity correlations: theory and practice,” ACM Transactions on Graphics, vol. 40, no. 3, pp. 1–22, 2021.
- 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,” Optics Express, vol. 23, no. 10, pp. 505–516, 2015.
- M. Hofer, C. Soeller, S. Brasselet, and J. Bertolotti, “Wide field fluorescence epi-microscopy behind a scattering medium enabled by speckle correlations,” Optics Express, vol. 26, no. 8, pp. 9866–9881, 2018.
- E. Betzig, G. H. Patterson, R. Sougrat et al., “Imaging intracellular fluorescent proteins at nanometer resolution,” Science, vol. 313, no. 5793, pp. 1642–1645, 2006.
- K. F. Tehrani, J. Xu, Y. Zhang, P. Shen, and P. Kner, “Adaptive optics stochastic optical reconstruction microscopy (AO-STORM) using a genetic algorithm,” Optics Express, vol. 23, no. 10, pp. 677–692, 2015.
- D. Wang, S. K. Sahoo, X. Zhu, G. Adamo, and C. Dang, “Non-invasive super-resolution imaging through dynamic scattering media,” Nature Communications, vol. 12, no. 1, p. 3150, 2021.
- K. B. Petersen and M. S. Pedersen, “The Matrix Cookbook,” Technical University of Denmark, vol. 7, no. 15, p. 510, 2008.
- C. Dong, C. C. Loy, K. He, and X. Tang, “Image super-resolution using deep convolutional networks,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 38, no. 2, pp. 295–307, 2016.
- S. Ruder, “An overview of gradient descent optimization algorithms,” http://arxiv.org/abs/1609.04747.
- R. G. Baraniuk, “Compressive sensing,” IEEE Signal Processing Magazine, vol. 24, no. 4, pp. 118–121, 2007.
- R. Horisaki, Y. Ogura, M. Aino, and J. Tanida, “Single-shot phase imaging with a coded aperture,” Optics Letters, vol. 39, no. 22, pp. 6466–6469, 2014.
- Y. Nishizaki, R. Horisaki, K. Kitaguchi, M. Saito, and J. Tanida, “Analysis of non-iterative phase retrieval based on machine learning,” Optical Review, vol. 27, no. 1, pp. 136–141, 2020.
- C. A. Metzler, F. Heide, P. Rangarajan et al., “Deep-inverse correlography: towards real-time high-resolution non-line-of-sight imaging,” Optica, vol. 7, no. 1, pp. 63–71, 2020.
Copyright © 2022 Yuto Endo et al. Exclusive Licensee Zhejiang Lab, China. Distributed under a Creative Commons Attribution License (CC BY 4.0).