1. Introduction

Mapping the three-dimensional (3D) morphology of an object has been an appealing subject for a wide range of disciplines including engineering, biomedicine, and life science. At a microscopic scale, several technologies such as optical coherence tomography (OCT) [1–3], confocal microscopy [4–6], multiphoton imaging [7–9], and digital holographic techniques [10–12] have proven to be useful, especially for observing the 3D structures of biological specimens. For macroscopic objects, focus has mainly been given to the fast reconstruction of real-world objects for various practical applications. Examples include the recognition of motion [13, 14] and navigation of automated systems [15, 16]. One approach for 3D macroscale imaging involves determining the time-of-flight of light returning from target objects by measuring the phase delay of the amplitude-modulated light [17–20]. This method is advantageous for fast measurements because a full field image can be acquired over multiple depths in a single recording. However, it is unsuitable for precise depth discrimination due to the limited bandwidth of amplitude modulation. The other approach is to use spatially modulated illumination for observing 3D shapes of macroscopic objects [21, 22]. Although this method allows for high spatial resolution relative to the time-of-flight method, the associated depth resolution can still be low because the spatial frequency of the illumination is limited by the low numerical aperture of the imaging optics.

In order to obtain both high spatial resolution and depth-selectivity, interferometric methods based on the short coherence length of light sources have been reported. Examples include white-light OCT [23, 24], extended focus method [25], structured illumination [26, 27], and phase-shifting or off-axis digital holography [28–30]. While these methods feature superior depth discrimination, they have mostly been applied to microscale objects. In addition, when objects are located behind a scattering medium, the quality of image reconstruction is significantly degraded due to the absence of the multiple scattering rejection mechanism. To resolve the multiple scattering problem in the macroscale imaging, several interesting methods have been introduced in which multiple light scattering was used for the image reconstruction. An image of a large-scale object was reconstructed by measuring the time-of-flight of the reflected light from a diffusive wall [31]. A wavefront shaping method has also been used to correct images obscured by scattering layers [32, 33] or control light through scattering media [34–36]. However, the method in [31] showed low spatial resolution due to the under-determined nature of image reconstruction and the method in [32] requires pre-calibration of the scattering layer. Imaging a human-sized object through a flame was reported using the multiple acquisitions of holographic images [37]. Similarly, we demonstrated high depth-resolution macroscopy capable of mapping the 3D morphology of an object hidden behind a scattering medium by using time-gated digital holographic imaging [38]. However, the enhancement of the image contrast was capped due to the incoherent averaging of multiple images.

For the effective suppression of multiple light scattering, we recently developed a unique method named collective accumulation of single-scattered waves (CASS) microscopy, in which both time-gated detection and aperture synthesis were used to collectively accumulate single-scattered waves [39]. The method was applied to imaging microscale objects embedded deep within a scattering media such as biological tissues. In this letter, we present CASS microscopy extended to the high depth-resolution imaging of a macroscopic object hidden behind a scattering layer. We first acquired depth-selective images for a macroscopic object using time-gated detection by interferometric measurement and then accumulated the images taken at different angles of illumination such that the single-scattered waves, i.e. the waves scattered only once by the object, were added in-phase. This combination resulted in the improvement of spatial resolution as well as the enhancement of signal to multiple scattering background ratio for the depth-selective imaging. We demonstrated 3D mapping of a macroscopic object through a 13-scattering-mean-free-path thick scattering layer, where conventional digital holographic imaging failed to work, with the depth resolution of 400 μm and view field of 30$\times$30 mm².

2. Methods and experimental setup

The schematic of the imaging system is presented in Fig. 1. The system is based on Mach-Zehnder interferometry using a diode laser (Thorlabs, LP637-SF70) as a light source. The center wavelength of the laser output was λ_c = 637 nm and its coherence length was 400 μm, which was significantly smaller than the thickness of the macroscopic object used in the experiment. For the interferometric detection, the output beam from the laser was split into a sample and a reference beam at the beam splitter (BS1). A two-axis galvanometer mirror was placed in the sample beam path to steer the incident angle of the sample beam to target objects. The object images were acquired at 200 different incidence angles covering a numerical aperture (NA) of ${\alpha }_{ill}=0.002$ at the sample plane. A reflection signal from the object was collected by L2 and directed to a camera (Lumenera, LM135, 1392$\times$1040 pixels, 4.65 μm pixel pitch) by a plate beam splitter BS2, which was placed between L2 and L3. The magnification from the sample plane to the camera was 1/5. The spatial-frequency bandwidth of the object spectrum collected by L2 was adjusted using the iris located at the focal plane of the lens. Specifically, the effective numerical aperture of L2 was set to ${\alpha }_{col}=0.0025$, which corresponds to the spatial resolution of 260 μm for the conventional coherent imaging. In the reference beam path, two mirrors were placed in a retro-reflection geometry on a scanning mirror for the fine-tuning of the path length. The scanning step of the reference arm for 3-D imaging was 200 μm. The reference beam was combined with the sample beam at the beam splitter (BS3) with an off-axis configuration. The interference image was then measured by the camera. The recorded interferograms were processed into complex field images providing phase and amplitude information using the standard Hilbert transform method [40].

 

Fig. 1 Experimental setup. Laser diode: 637 nm in wavelength. BS1, 3: cube beam splitters, L1 – L5: lenses, BS2: plate beam splitter, CCD: Camera, scattering layer: PDMS mixed with ZnO particles. The thickness of the scattering layer was 720.0 ± 0.5 $\text{μ}$m, and its scattering mean free path was $l_s$ = 55.9 ± 4.6 $\text{μ}$m.

Download Full Size | PDF

For a macroscopic scale test target, a Julian plaster model was used in the experiment. The height and width of the model were 6 cm and 3 cm, respectively, and its thickness along the axial direction was about 2 cm. In the experiment, 3$\times$3 cm² area near the forehead of the model object was imaged. The size of the view field was limited by the number of detector pixels. As a scattering layer, a mixture of Polydimethylsiloxane (PDMS) and Zinc-Oxide (ZnO) particles was used. The scattering mean free path l_s was determined by a separate experiment in which the transmittance of the ballistic photons was measured as a function of the thickness of the scattering layer. Exponential curve fitting to the measured transmittance led to the finding of the decay constant, which corresponds to the l_s of 55.9 ± 4.6 μm. We also determined the transport mean free path l_t as 383.0 ± 0.4 μm by measuring the total transmittance of light as a function of the thickness. The thickness of the layer for imaging was prepared as 720.0 ± 0.5 μm, which corresponds to approximately 13 times the l_s. This scattering layer was placed between the test target and L2 to cover the entire field of view.

3. Experimental setup and results

Representative single-shot images acquired for the Julian plaster model are presented in Fig. 2. Figure 2(a) shows a single raw interference image taken without the scattering layer. Although the entire shape of the object could be seen, the interference pattern was only found within a very limited region due to the short coherence length of the laser. The area indicated by a red box is magnified in the inset in Fig. 2(a) to clearly visualize the interference fringes. After the application of the Hilbert transform, a band-shaped amplitude distribution was observed as shown in Fig. 2(b). This shows the contour of the target object having the same depth as that set by the position of the reference mirror. Figure 2(c) shows a single raw interference image for the object taken after the insertion of the scattering layer. In contrast to the case without the scattering layer, the shape of the object could not be resolved in the raw image due to the multiple scattering generated by the scattering layer. In addition, the interference fringe patterns became extremely weak despite matching the optical path length of the reference wave. Fringe patterns were almost invisible in the magnified image in the inset in Fig. 2(c). Consequently, the contour of the object was not visible in the amplitude map of the Hilbert-transformed image in Fig. 2(d).

 

Fig. 2 Single-shot images taken with and without a scattering layer. (a) Raw interference image of the Julian plaster model in the absence of a scattering layer. Sample beam was normally incident to the target object. Inset: magnified image of the content in the small red box. (b) Reconstructed image using interference pattern of the single-shot raw image in (a). (c) and (d): Same as (a) and (b), respectively, but with the insertion of a scattering layer. Scale bar, 5 mm. Color bars in (a) and (c) indicate intensity in arbitrary units and those in (b) and (d) indicate amplitude in arbitrary units.

Download Full Size | PDF

Since individual single-shot images could not resolve object structures due to multiple scattering, we employed coherent aperture synthesis to suppress the effect of multiple light scattering. We acquired multiple images for various illumination angles by the scanning of the galvanometer mirror (GM) in the sample beam path. The scanning was performed so that 200 illumination angles uniformly covered the numerical aperture of ${\alpha }_{ill}$. First, the amplitude maps in the spatial frequency and real spaces are shown in Figs. 3(a) and 3(d), respectively, for the case of normal illumination. As shown in Fig. 3(d), it is almost impossible to detect the object features in a single-shot image. Due to the normal illumination, the detected spectrum in Fig. 3(a) was distributed around the center in k-space. The radius of the spectrum denoted by the black dotted circle was $k_0{\alpha }_{col}$. Due to the multiple-scattered waves induced by the scattering layer, the detection NA was fully filled with random speckles.

 

Fig. 3 Aperture synthesis of multiple images taken at a fixed depth after the insertion of a scattering layer. (a) and (b) Amplitude maps of reflected waves in the $\left(k_x,k_y\right)$ space under a normal illumination and an oblique illumination corresponding to incident wave vector of ${\stackrel{\rightharpoonup }{k}}^i$ = $(k_0{\alpha }_{ill}$, 0), respectively. Scale bar, $k_0{\alpha }_{col}$. The radius of the black and blue circles is $k_0{\alpha }_{col}$. (d) and (e) Amplitude maps of the 2D inverse Fourier transformed images of (a) and (b), respectively. Scale bar, 5 mm. Color bar, amplitude in arbitrary units. (c) Angular spectrum after the aperture synthesis of all 200 images taken at various illumination angles. The radius of the red circle is $1.8\times k_0{\alpha }_{col}$. (f) 2D inverse Fourier transformed image of (c), showing the synthetic aperture image reconstructed by adding all 200 complex images. All angular spectra are represented in log scale for better visibility.

Download Full Size | PDF

Figures 3(b) and 3(e) show the amplitude maps of the acquired image in the spatial frequency and real spaces, respectively, for the spatial frequency of illumination ${\overrightarrow{k}}^i=\left(k_0{\alpha }_{ill},0\right)$. Similar to the normal illumination shown in Fig. 3(a), multiple-scattered waves filled the entire collection bandwidth indicated by the blue dotted circle in Fig. 3(b), and the object features could not be resolved.

We synthesized all the acquired images based on the object spectrum. In other words, the angular spectrum of each image was shifted by the wave vector of the illumination. For instance, as shown in Fig. 3(c), the i-th angular spectrum in Fig. 3(b) was added to that in Fig. 3(a), after the shift of the spectrum by $-k^i$. This is the amount of the i-th oblique illumination which was obtained by the pre-calibration for the illumination angle. During this aperture synthesis process, the single-scattered waves from the object were added in phase because the object spectrum was made static for different illumination angles. In contrast, the multiple-scattered waves have random phases and have no correlation among them for different illumination angles. Thus the multiple-scattered waves were added with almost no cross terms, resulting in the incoherent addition. In addition, the spectral bandwidth for imaging was enlarged by the amount of angular coverage of the GM scanning. The resulting bandwidth of the synthesized image was $k_0\left({\alpha }_{col}+{\alpha }_{ill}\right)$, which is the radius of the red circle in Fig. 3(c). Figure 3(f) shows the final image obtained by the 2D inverse Fourier transform shown in Fig. 3(c). Due to the coherent accumulation of the single-scattered waves and incoherent addition of multiple-scattered waves, the contour line of object features associated with the coherence length of the laser could be clearly resolved, unlike the single-shot image.

To demonstrate the effectiveness of the signal enhancement of the time-gated aperture synthesis, we also reconstructed the object image by the incoherent addition. All the 200 acquired complex images were first converted into intensity images and then added with no consideration of their phase. We compared this reconstruction result with that produced by the aperture synthesis. Figures 4(a) and 4(b) show the intensity maps for the object reconstructed at the same depth by the incoherent addition and the aperture synthesis, respectively, without the scattering layer. As expected, the signal strength was significantly enhanced and the object features were clearly visible in both images. Although the two images appear to be similar, they significantly differ in terms of the signal enhancement. Compared to the single-shot image in Fig. 2(b), the signal growth factors were about 170 and 21,000 for the incoherent addition and aperture synthesis, respectively. This is because the signal increases linearly for the incoherent addition, while it undergoes quadratic growth for the case of the aperture synthesis. Slight mismatches in the enhancement factors might be due to the experimentally uncontrollable conditions such as the angular variation of the object signal and the system vibration.

 

Fig. 4 Comparison between incoherent imaging and synthetic aperture imaging. (a) and (b) Incoherent image and synthetic aperture image, respectively, without the scattering layer. Scale bar, 5 mm. Color bars indicate intensity in arbitrary units, but in the same scale for (a) and (b). (c) and (d) Same as (a) and (b), respectively, but with the insertion of the scattering layer. Color bars, intensity in arbitrary units.

Download Full Size | PDF

We then examined the case in which the scattering layer was positioned. Figures 4(c) and 4(d) present the reconstructed images by the incoherent addition and the aperture synthesis, respectively. Unlike the case without the scattering layer where both methods can produce clear object images, only the aperture synthesis in Fig. 4(d) could resolve the target structure. The reduction of the image contrast with the scattering layer compared to the case of no multiple scattering can be explained by the dependence of signal growth mechanism to the ratio of single- to multiple-scattered waves. The detailed explanation is given in the following.

Statistically, the signal intensity S (the average intensity of all images for the region where the object signal exists) of the recorded single-shot images can be written as S = s₀ + m₀, where s₀ is the intensity of single-scattered waves containing the correct object information and m₀ is the intensity of multiple-scattered waves which are mostly noise. In the incoherent addition, both s₀ and m₀ increase linearly, proportional to N, the number of taken images. After the addition, the intensity for the reconstructed image can be written as S_inc = Ns₀ + Nm₀. Thus, the ratio for the signal intensity of the reconstruction to that of the single-shot image is S_inc/S = N. This is the same as that for the case without multiple scattering.

We assessed this point in the measured images. The region of interest where the object existed was selected from an appropriate threshold in the single-shot image taken without the scattering layer, shown in Fig. 2(b). The average of the intensity within this region was taken as S for the corresponding image. We used the same region for obtaining the signal intensity S for other images taken at different illumination angles. As shown in Figs. 5(a) and 5(b), S_inc grew almost linearly in proportion to N, irrespective of the presence of the scattering layer. The final growth factors were about 170 and 200 with and without the scattering layer, respectively, when N = 200. By fitting the growth curves in log-log scale, the slopes were determined as 0.961 ± 0.001 and 0.992 ± 0.001 for the cases of with and without the scattering layer, respectively. These values are fairly close to 1, which confirms the linear growth of the signal intensities by the incoherent addition process. Although the total intensity increases by N, the ratio of the single-scattered component to the multiple-scattered component becomes γ_inc = Ns₀/Nm₀ = s₀/m₀ = γ₀, where γ₀ is the ratio of single- to multiple-scattered components of individual single-shot images. Therefore, the object signal relative to the background remains statistically unchanged after the reconstruction by the incoherent addition. Since the object feature was unresolvable in the single-shot images shown in Fig. 2(d), the multiple scattering is expected to dominate the single scattering, i.e. γ₀≪1 in the presence of the scattering layer. This explains why the target object could not be revealed by the incoherent addition in Fig. 4(c).

 

Fig. 5 Signal growth versus the number of incidence angles used for image reconstruction. S_syn and S_inc were normalized by the average signal intensity S. The red curve represents the results obtained by the aperture synthesis, and the blue curve represents those by the incoherent addition for the cases without (a) and with (b) the scattering layer. Dashed and dashed-dot lines show growth by N² and N for the visual guidance for the ideal aperture synthesis and incoherent addition, respectively. Graphs are represented in log-log scale.

Download Full Size | PDF

Unlike the incoherent addition, the single-scattered waves increases by $N^2$ in the aperture synthesis due to the coherent superposition of the signal fields. On the other hand, the multiple scattering background still grows by N because of no phase correlation among the components of multiple-scattered waves. As a result, the total intensity of the reconstructed image at the target object is written as S_syn = N²s₀ + Nm₀ in the aperture synthesis and the signal growth factor S_syn/S becomes (N²s₀ + Nm₀)/(s₀ + m₀). When multiple scattering does not occur, this factor shows quadratic growth with respect to N. This was confirmed in Fig. 5(a) where the signal grew up to 21,000 when N = 200. This signal growth curve was close to that of the guideline (black dashed line), representing N² growth. The slope of signal growth in the log-log plot was fitted to be 1.845 ± 0.004, which is close to the theoretically expected value of 2. The discrepancy was found to increase as N increases. This is due to the phase fluctuations induced by the system vibration, which causes imperfect addition of measured complex fields. The reduction of the measured signal at high illumination angles was also responsible for the decrease of the signal growth rate for large N.

In the presence of a scattering layer, the growth of the signal at the target is lower than that in the ideal case. In the extreme case when the multiple scattering dominates single scattering, i.e. γ₀ ≪ 1, signal growth can be approximated to S_syn/S ≈N²s₀/m₀ + N = N²γ₀ + N. This can be further approximated as S_syn/S ≈ N²γ₀ when the number of measurements is large enough such that N≫1/γ₀. The signal growth in turn is smaller than N² by a factor of γ₀. Even with the additional factor, the signal grows quadratically by N², similar to the case of no scattering. To see the N²-dependence, we fitted S_sys in Fig. 5(b) in the range of large N (N>30), where N²s₀ is sufficiently larger than Nm₀. The slope of the signal growth in the log-log plot was fitted to be 1.724 ± 0.003, which is close to the N² growth. The slope was smaller than the case with no scattering layer, and this can be attributed to the poor SNR of the single-shot images, which may affect the precision of the aperture synthesis.

In the experimental data shown in Fig. 5(b), S_syn/S reaches about 4,100 after the coherent aperture synthesis with N = 200. This implies that γ₀ was about 0.1 under our experimental conditions, confirming why the object was invisible in each single-shot image. On the contrary, the ratio of the single- to multiple-scattered components was increased to γ_syn = N²s₀/Nm₀ = Nγ₀ = 20. After the aperture synthesis, the single scattering signal became considerably larger than the multiple scattering signal, i.e. ${\gamma }_{syn}\gg 1$, so that the object features could be clearly resolved with the enhanced signal contrast.

We repeated the same measurements for different depths, and the 2D contour image for each depth was combined to form point clouds of the target objects in Fig. 6. The imaging depth was chosen by translating the position of the reference mirror. Due to the restriction on the depth of focus of the imaging lens, the depth scanning range was limited to 8.8 mm from the forehead, covering part of the volume of the plaster model.

 

Fig. 6 Comparison of 3D reconstructed images between incoherent imaging and synthetic aperture imaging. (a) and (b) 3D reconstruction images using incoherent imaging and synthetic aperture imaging, respectively, in the absence of the scattering layer. (c) and (d) Same as (a) and (b), respectively, but in the presence of the scattering layer.

Download Full Size | PDF

Figures 6(a) and 6(b) show the 3D maps generated by the incoherent addition and the aperture synthesis, respectively, when no scattering layer was present. Similar to the cases in Figs. 4(a) and 4(b), the object shape could be clearly visualized in 3D for both cases. For the case when the scattering layer was present, 3D maps produced from the incoherent addition and the aperture synthesis are presented in Figs. 6(c) and 6(d), respectively. Unlike the case in Fig. 6(a), the incoherent addition could not show any object features due to the poor signal-to-noise ratio (SNR) caused by the multiple scattering. The aperture synthesis, however, successfully generated noticeable object distribution, as shown in Fig. 6(d). Although the image contrast was somewhat reduced due to the presence of the multiple scattering, the aperture synthesis successfully revealed the object structure.

4. Conclusion

In conclusion, we developed a high-resolution imaging method which is capable of retrieving the 3D morphology of a macroscopic object hidden behind a scattering layer. Multiple time-gated reflection images were taken at various angles of illumination and the object information was enhanced by coherently synthesizing the object apertures. This resulted in the collective accumulation of single-scattered waves from the target object and consequently enhanced the SNR beyond that achieved by the method of incoherent addition. With the capabilities of enhancing the single-scattered waves, higher spatial resolution, and superior imaging sensitivity, our method can potentially be used for precision 3D measurements in a demanding macroscopic world where objects are obscured by scattering media.