Single-shot non-invasive three-dimensional imaging through scattering media

TL;DR

This paper introduces a novel single-shot 3D imaging technique through scattering media using a 3D memory effect, enabling 3D object recovery from a single speckle image without multiple exposures.

Contribution

It presents a new computational method combining speckle correlation and phase retrieval for single-shot 3D imaging through scattering media, improving upon multi-shot approaches.

Findings

Successfully demonstrated lensless 3D imaging through scattering media.

Achieved 3D object recovery from a single speckle image.

Compared favorably with previous multi-shot methods.

Abstract

We present a method for single-shot three-dimensional imaging through scattering media with a three-dimensional memory effect. In the proposed computational process, a captured speckle image is two-dimensionally correlated with different scales, and the object is three-dimensionally recovered with three-dimensional phase retrieval. Our method was experimentally demonstrated with a lensless setup and was compared with a multi-shot approach used in our previous work [Y. Okamoto, et al., Opt. Lett. 44, 2526-2529 (2019)].

Equations

$$i({\mathbf{r}}_i)=\int h({\mathbf{r}}_i-{\mathbf{r}}_o)o({\mathbf{r}}_o){\text{d}}^3{r}_o,$$

$$s=(z_o+z_i)/z_o.$$

$$i{\star }_{\text{3D}}i$$

$$=(h{\star }_{\text{3D}}h){\otimes }_{\text{3D}}(o{\star }_{\text{3D}}o)$$

$$\approx o{\star }_{\text{3D}}o$$

$$={\mathcal{F}}_{\text{3D}}^{-1}[|O{|}^2],$$

$${c^{k}=}\left\{\begin{array}[]{ll}i_{\text{2D}}^{0}\star_{\text{2D}}i_{\text{2D}}^{k},&\text{for}\leavevmode\nobreak\ \leavevmode\nobreak\ k\geq 0\\ i_{\text{2D}}^{-k}\star_{\text{2D}}i_{\text{2D}}^{0},&\text{for}\leavevmode\nobreak\ \leavevmode\nobreak\ k<0\end{array}\right.$$

$$o_{n}(\bm{r}_{o})=\left\{\begin{array}[]{ll}o^{\prime}_{n}(\bm{r}_{o}),&\text{for}\leavevmode\nobreak\ \leavevmode\nobreak\ \bm{r}_{o}\notin\Gamma\\ 0,&\text{for}\leavevmode\nobreak\ \leavevmode\nobreak\ \bm{r}_{o}\in\Gamma\end{array}\right.$$

$$o_{n}(\bm{r}_{o})=\left\{\begin{array}[]{ll}o^{\prime}_{n}(\bm{r}_{o}),&\text{for}\leavevmode\nobreak\ \leavevmode\nobreak\ \bm{r}_{o}\notin\Gamma\\ o_{n-1}(\bm{r}_{o})-\beta o^{\prime}_{n}(\bm{r}_{o}),&\text{for}\leavevmode\nobreak\ \leavevmode\nobreak\ \bm{r}_{o}\in\Gamma\end{array}\right.$$

$$s_{\text{com}}^m=\frac{(z_o-m{\Delta }_{z_o}+z_i)/(z_o-m{\Delta }_{z_o})}{(z_o+z_i)/z_o},$$

$$|s_{\text{com}}^{M-1}-1|d\le \delta ,$$

Full text

††journal: ao

1 Introduction

Imaging through scattering media has been long studied for biomedical imaging, astronomical imaging, and so on [1, 2, 3, 4]. Recently, optical sensing and control through strongly scattering media, which are difficult to handle in conventional approaches relying on the existence of non-scattered light, have attracted interest in the field of optics and photonics. The rapidly growing computational power and the improved performance of optical elements for light control drive this area, and various methods have been reported. These methods are categorized into three types: feedback-based, inversion-based, and correlation-based [5, 6].

The feedback-based approach utilizes wavefront shaping behind or inside scattering media with an iterative feedback process based on an optimization algorithm [7, 8, 9, 10, 11, 12, 13, 14, 15]. Issues with the feedback-based approach are the large number of feedback loops and the need for a probing process to measure the focusing state in the every loop. The inversion-based approach, including time reversal and phase conjugation, senses and controls the optical distribution through scattering media by taking the inverse of a transmission matrix expressing the scattering process [16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30]. These methods realize single-shot imaging and focusing without any feedback process. However, they need to probe the whole or part of the transmission matrix before the imaging and focusing stage.

The correlation-based approach exploits the shift invariance of speckles, which is called the memory effect [31, 32, 33, 34, 35, 36, 37, 38, 39]. In the correlation-based methods, an autocorrelation process is used for removing speckles and exposing object signals in captured images. An advantage of the correlation-based approach is the lack of a need for the probing process, which is a drawback of the previous two approaches. As a result, the correlation-based approach has realized non-invasive imaging through scattering media. This approach has been recently extended to a three-dimensional case [40, 41, 42, 43].

Drawbacks with the methods for correlation-based three-dimensional imaging through scattering media include the need for capturing multiple images and/or invasive optical processes. Thus, they are difficult to apply to imaging of dynamical and practical scenes. Here we present a method for tomographic reconstruction of a three-dimensional object from a single speckle image captured through scattering media without any invasive or probing process. We demonstrated the method experimentally with a lensless setup. Our method enhances the possibility and practicability of imaging though scattering media for a wide range of applications, including biomedicine, security, and industry.

2 Method

The optical setup in the proposed method is shown in Fig. 1. A three-dimensional object $o$ is illuminated with spatially incoherent illumination through a first diffuser, and the light scattered by a second diffuser is captured by a lensless image sensor as $i$. The relationship between the object and the captured speckle is three-dimensionally shift-invariant, based on the three-dimensional memory effect [41]. Then, the imaging process is written with a scattering impulse response $h$ as

[TABLE]

where $\bm{r}_{i}=(x_{i},y_{i},z_{i})$ are the spatial coordinates in the object space, and $\bm{r}_{o}=(x_{o},y_{o},z_{o})$ are the spatial coordinates in the sensor space, respectively. The $x$- and $y$-axes are lateral to the image sensor, and the $z$-axis is axial to the image sensor. The origin of these coordinates is the center of the second diffuser. The impulse response $h$ is laterally random and axially scaled with a scaling factor $s$, which is written as

[TABLE]

In our previous work on multi-shot three-dimensional imaging through scattering media [43], the image sensor is axially scanned to observe the three-dimensional speckle distribution $i$, which is three-dimensionally autocorrelated to remove the impact of the scattering process $h$ as

[TABLE]

where $\star_{\text{3D}}$ is the three-dimensional correlation, $\otimes_{\text{3D}}$ is the three-dimensional convolution, $\mathcal{F}_{\text{3D}}$ denotes the three-dimensional Fourier transform, and $O$ is the three-dimensional transform of $o$. Here, $h\star_{\text{3D}}h$ is approximated by the delta function due to the laterally random and axially scaled distribution of $h$ [33, 34]. The object signal $o$ is recovered by a phase retrieval process for $|O|^{2}$ [44].

In this study, as shown in Fig. 2, a three-dimensional object is tomographically reconstructed from a single captured speckle image $i_{\text{2D}}$, which contains differently scaled impulse responses $h$ depending on the object distances $z_{o}$, using the optical setup in Fig. 1. The captured speckle image is computationally and laterally scaled multiple times, instead of multiple measurements made while axially scanning the image sensor as in the previous work [43]. This computational scaling process mimics the optical scaling process in Eq. (2) by axial scanning, where the scaling factor is monotonically increased or decreased. This process is exploited for searching relative scales of the impulse responses on the original speckle image through the following correlation process, where high correlations appear if the original and scaled impulse responses are coincident, and vice versa.

The scaled speckle images are computationally and laterally correlated as follows:

[TABLE]

where $\star_{\text{2D}}$ is the two-dimensional correlation, the superscript $k=-M+1,-M+2,\cdots,M-2$ denotes the index for the lateral correlations, $M$ is the number of scaled images, $c^{k}$ is the $k$-th lateral correlation result, $i_{\text{2D}}^{m}$ is the $m$-th scaled speckle image ($i_{\text{2D}}^{0}=i_{\text{2D}}$), and the superscript $m=0,1,\cdots,M-1$ denotes the index of the scaling factor. In this study, $o\star_{\text{3D}}o$ in Eq. (5) is approximated by $c^{k}$. The object $o$ is three-dimensionally reconstructed from the correlation result $c^{k}$ based on three-dimensional phase retrieval.

The phase retrieval process in this study follows the previous work on speckle-correlation-based imaging through scattering media [33, 34, 43]. Here two algorithms, which are called the error reduction algorithm and the hybrid input-output algorithm [44], are sequentially performed, as shown in Fig. 2. The second algorithm is a modified version of the first one, and they use the loop in the iterative processes shown in Fig. 3. The difference between them is the process of constraints. First, the common aspects are described, and then the differences are explained.

The iterative process of the error reduction algorithm and the hybrid input-output algorithm is shown in Fig. 3 and is described as follows:

The object’s Fourier spectrum is initially given with a three-dimensional random phase $\theta_{n}$ by $O_{n}=(|O|^{2})^{(1/2)}\exp(j\theta_{n})$, where $j$ is the imaginary unit, $|O|^{2}$ is the three-dimensional Fourier transform of $c^{k}$ in Eq. (7), and the subscript $n$ is the counter of the iteration, which is set to one. 2. 2.

$O_{n}$ is three-dimensionally inverse Fourier transformed, and the result is set as the intermediately estimated object $o^{\prime}_{n}$. 3. 3.

$o^{\prime}_{n}$ is rectified with some constraints, which are described in the following paragraphs, and the result is the estimated object $o_{n}$ at the $n$-th iteration. 4. 4.

$o_{n}$ is three-dimensionally Fourier transformed, and the result is set as the intermediately estimated object’s Fourier spectrum $O^{\prime}_{n}$. 5. 5.

The argument of $O^{\prime}_{n}$ is extracted, and it is used as a replacement of the argument $\theta_{n}$ of the estimated object’s Fourier spectrum $O_{n}$, where the counter $n$ is incremented by one.

Steps 2–5 are iterated.

The rectifying process with constraints at Step 3 is the difference between the error reduction algorithm and the hybrid input-output algorithm [44]. In the error reduction algorithm, the intermediately estimated object $o^{\prime}_{n}$ at the $n$-th iteration is updated with the following rule:

[TABLE]

where $\Gamma$ is the set of all spatial positions $\bm{r}_{o}$ which violate the constraints. In the hybrid input-output algorithm, the updating rule is as follows:

[TABLE]

where $\beta$ is a feedback parameter.

In the reconstruction according to the proposed scheme, as shown in Fig. 2, the hybrid input-output algorithm is performed first. In this algorithm, the feedback parameter $\beta$ is decreased from 2.0 to 0.0 in intervals of 0.05, and the loop of Fig. 3 is iterated ten times for each $\beta$. Then, the error reduction algorithm is performed by using the result of the hybrid input-output algorithm as the initial estimate at Step 1, and the loop of Fig. 3 is iterated five hundred times. The constraints used here are realness, non-negativity, and the range of pixel intensities. Realness and non-negativity are introduced by using spatially incoherent illumination, such as a light emitting diode (LED) or fluorescence [33, 34]. The range of pixel intensities suppresses some reconstruction artifacts [43]. This phase retrieval has trivial ambiguities of the spatial shift and the conjugate inversion, which have been studied in the literature on phase retrieval [45].

3 Experimental demonstration

The proposed method was demonstrated with three-dimensionally arranged point sources fabricated by a 3D printer (M2030TP manufactured by L-DEVO). The object had three levels, where the step height was 0.5 cm, and holes with a diameter of 0.4 mm were located at different positions on each level, as shown in Fig. 4. One hole was made on the front level, two diagonally arranged holes were made on the middle level, and two vertically arranged holes were made on the back level, respectively. As shown in Fig. 1, the object was located between two diffusers and it was illuminated with an incoherent LED (M565L3 manufactured by Thorlabs, nominal wavelength: 565 nm, full width at half maximum of spectrum: 103 nm) through a bandpass filter (578NM X 16NM 25MM manufactured by Edmund Optics, central wavelength 578 nm, full width at half maximum of spectrum: 22 nm) and a diffuser (LSD5PC10-5 manufactured by Luminit). The distance between the first diffuser and the object was 75 mm. Light passing through object was scattered by another diffuser (LSD20PC10-5 manufactured by Luminit) and captured by a monochrome image sensor (hr29050MFLGEA manufactured by SVS-Vistek, pixel count: $4384\times 6576$, pixel pitch: $5.5\times 5.5$ µm) without any imaging optics. The distance ($z_{o}$) between the object and the second diffuser was 92 mm, and the distance ($z_{i}$) between the second diffuser and the image sensor was 25 mm, as shown in Fig. 1.

The captured speckle image is shown in Fig. 4, where the central $600\times 600$ pixel area was clipped for visualization purposes. After background compensation, the captured image was scaled with a computational scaling factor $s_{\text{com}}^{m}$ based on Eq. (2) as follows:

[TABLE]

where $\Delta_{z_{o}}$ is the axial resolution of the object space, which was set to 0.5 mm, corresponding to the step height of the object. The number of scaled speckle images ($M$) was set to six. The scaled speckle images were laterally correlated as in Eq. (7). The central $800\times 800$ pixel areas of the correlations were clipped and laterally down-sampled by a factor of four to reduce the computational cost in the next phase retrieval process. Biases of the correlations were equalized with the average values outside the central areas. The result of this correlation process is shown in Fig. 4, where the central $150\times 150$ pixel areas were clipped, and contrast enhancement was applied for visualization purposes.

The three-dimensional phase retrieval process was applied to the correlations, and the results are shown in Fig. 4, where the ambiguities of the spatial shift and the conjugate inversion were manually compensated, and the central $150\times 150$ pixel areas were clipped for visualization purposes. The holes were three-dimensionally recovered. The reconstruction result of the multi-shot approach with axial scanning of the image sensor in our previous work [43] is shown in Fig. 4. These results were comparable, although some interference from other planes in the single-shot approach was stronger than that in the multi-shot approach.

An issue with the single-shot approach compared with the multi-shot one is distortion of the autocorrelation through the scaling process. The distortion should be smaller than the lateral resolution of speckle correlations. The limitation of the scaling factor in Eq. (10) is calculated as

[TABLE]

where $d$ is the object size on the image sensor and is estimated as half of the correlative area, and $\delta$ is the resolution of speckle correlations and is $\sqrt{2}$-times larger than the grain size of the speckles, assuming a Gaussian distribution. In the case of the experiment, $d=300\leavevmode\nobreak\ \text{pixels}$ and $\delta=\sqrt{2}\times 17\leavevmode\nobreak\ \text{pixels}$ [43], so the above limitation was satisfied.

4 Conclusion

We proposed single-shot three-dimensional imaging through scattering media based on speckle correlation. An object is three-dimensionally reconstructed from a single speckle image with a scaling process, a correlation process, and a phase retrieval process. The proposed method was experimentally demonstrated with three-dimensionally arranged point sources between diffusers. The result was comparable to the multi-shot case reported previously.

Our method simplifies the optical setup for three-dimensional imaging through scattering media. It is useful in a wide range of modalities, including lensless imaging, microscope imaging, and telescope imaging, and applications such as looking around corners and reflection-mode scattered imaging. Also, this method may provide interesting insights for three-dimensional imaging without the need to calibrate optical hardware or probe optical phenomena before an imaging stage.

Funding

This work was supported by JSPS KAKENHI Grant Numbers JP17H02799 and JP17K00233, JST PRESTO Grant Number JPMJPR17PB.