Scattering media widely exist in our daily life, which possess a non-uniform reflective and refractive index distribution. It disturbs the light rays coming from the imaging targets, which hinders the direct analysis of object information behind it using traditional optical systems. To overcome this challenge, nowadays, novel methods based on the memory effect (ME) [1–4] achieve non-invasive scattering imaging via speckle correlation [5–7]. Compared with other scattering imaging methods, such as the ballistic-light-based approach [8–13], wavefront shaping [14,15], and transmission matrix measurement [16–21], the speckle correlation technique realizes prior-free imaging only by traditional instruments and has the capability of quick imaging in currently inaccessible scenarios. However, the field-of-view (FOV) of this method is limited by the ME range. Prior-free imaging beyond the ME range is critical to seeing through scattering media.

To exceed the ME range, some techniques have been proposed by introducing the prior information of isolated imaging targets. Li et al. introduced the position prior of each imaging target during the point spread function (PSF) calibration [22]; Sahoo et al. used the wavelength prior of each imaging target during the PSF calibration [23]; and Guo et al. exploited the shape prior or the PSF closing to the imaging target to exceed the ME range [24]. Wang et al. proposed a dual-target non-invasive scattering imaging method via Fourier spectrum guessing and iterative energy constrained compensation [25]. However, this method can only be used for dual-target separation and the number of targets serves as the known prior information before reconstruction. Boniface et al. achieved non-invasive target localization beyond the ME by analyzing the speckle envelope of each target [26], but this localization method was not available when an unknown number of targets were illuminated simultaneously.

Here, we put forward a multi-target large FOV scattering imaging method based on the blind target position detection. It blindly detects the unknown number and positions of the isolated targets only using two multi-target near-field speckles captured at different imaging distances. The theoretical scaling relationship between two speckles is derived and demonstrated and the scaling centers correspond to the positions of imaging targets. Based on the theoretical derivation, scaling-vector-based target position detection is proposed, which recognizes the target position using the length and direction information of scaling vectors. After that, the autocorrelations can be calculated for speckle regions centered by the derived positions via the low-cross-talk region allocation strategy. Working with the modified phase retrieval algorithm [27] to select the optimal recovery with no prior information, especially for the autocorrelation with interference, the complete scene of multiple isolated targets exceeding the ME range is reconstructed. Experiments on a real imaging system demonstrate the effectiveness of the proposed algorithm in multi-target blind reconstruction through scattering media. Visually distinguishable reconstructions are experimentally achieved with the whole scene of multiple isolated targets exceeding the ME range. In principle, the exceeding times can be increased as long as the acquisition equipment can fully collect the scattered speckle in larger FOV. Further, we verified the accuracy for multi-target positioning, capability for multi-target blind reconstruction, and the experimental analysis about limitations and restraints of the proposed method in the end.

The principle of the method is depicted in Fig. 1. We excite multiple isolated targets, O1,O2,…,On, hidden behind a scattering medium, with a spatially incoherent illumination. Each imaging target falls within the ME range, while the spacing between any two targets is beyond. The large FOV speckle image at Distance 1, Id1, which consists of speckles produced by each target, Id11,Id12,…,Id1n, is captured by a 2D camera array. Certainly, the proposed imaging system is applied for the situation that the number and locations of the imaging targets cannot be seen directly from the captured speckle. This position information is hidden in Id1. To extract it, a second speckle Id2 is captured at another imaging distance, Distance 2. The ME range, Δθ, corresponds to the angle FOV of the diffuser, within which the points on the object plane produce random speckles with high correlation. For the spatially incoherent imaging system via speckle correlation, its ME range is constrained by Δθ≪λ/πL [5], where λ denotes the wavelength of the spatially incoherent light source and L is the effective thickness of the scattering media. Since the speckle patterns generated within the ME range are translational invariant, the large FOV speckle of multiple targets recorded by the camera array at Distance 1, Id1 in Fig. 1 can be formulated by (so as d2 for Id2) where * denotes the convolution operation; n denotes the total number of the imaging targets and Ok(k=1,2,…,n) are all the imaging targets; PSFd1k represents the translational invariant point spread function at Distance 1 corresponding to the point where Ok is located in the object plane; and Id1k denotes the speckle pattern generated by Ok with imaging distance d1. The translational invariance of the PSF can be further expressed by the envelope of each PSF, which is normally removed during the image processing and seen as the obstacle for reconstruction [5]. It varies with the distance away from the corresponding point of each PSF. Especially for multi-target scattering imaging, multiple envelopes couple the multi-target information with the hidden target position information, which increases the difficulty for imaging exceeding the ME. The intensity distribution of each PSF can be divided into two parts by where (x,y) is 2D coordinates in the sensor plane; (u,v) is 2D coordinates in the object plane, and (u=uk,v=vk) is the location of Ok in the object plane. C denotes the envelope property (or energy distribution of the PSF), which has high value if (x,y) is close to (uk,vk) with relatively low value if (x,y) is far away from (uk,vk) [1] and shows the same property for the PSF produced by a point-like source or the speckle generated by a small size object. Sd1k denotes the system response at Distance 1 corresponding to Ok after removing the envelope and the autocorrelation of Sd1k equals a sharp peak function [5]. Substituting Eq. (2) into Eq. (1), we have

From Eq. (3), the number, the shapes, and the locations of the imaging targets are coupled together in the captured speckle. This is the problem with an infinite number of solutions when reconstructing all the target shapes only via Id1, since there are so many unknown variables in Eq. (3). One effective way for simplifying this problem and reducing the solution space is to detect the total number of the objects n and the location of each object (uk,vk),k=1,2,…,n without any prior information. Thus, the blind target position detection algorithm is proposed and introduced in the following.

We exploit a second near-field speckle Id2 captured at another imaging distance for position detection. The theoretical scaling relationship between the two near-field speckles derived by us indicates that some areas (called scaling centers) in the near-field speckles do not scale with the imaging distance and these scaling-invariant areas correspond to the locations of imaging targets. To the best of our knowledge, no work has been proposed to investigate the relationship between two near-field speckles at different imaging distances. Inspired by the existing statistical analysis of the far-field speckles in 3D space [18,28,29], and the novel technique using speckle correlation to improve axial sectioning [30], without loss of generality, we first analyze the scaling relationship between two near-field on-axis PSFs. With an ideal pinhole at the optical axis in the object plane as the corresponding point where one imaging target (supposed as Ok) is located in the object plane, using the Fresnel diffraction formula, the field on the front surface of the scattering media, Us, is expressed by where we assume that the on-axis pinhole is located at (u=0,v=0); (xs,ys) is the 2D coordinates corresponding to the diffuser plane. Here, the scattering medium regarded as an unknown random 2D phase disturbance, TM(xs,ys), is introduced into the propagation model. Using the near-field Huygens diffraction theory, the light field on the first sensor plane, h1, can be expressed by where 1/f equals (1/d0+1/d1). The scattered PSF shown in Eq. (2) at Distance 1 plane, PSFd1k (so as d2 for PSFd2k at Distance 2 plane), is the square of the magnitude of h1 as

To describe the scaling relationship, we introduce the correlation function between on-axis PSFd1k and PSFd2k as where m is a constant value which represents that the spatial coordinates of PSFd2k are scaled to PSFd1k. In theory, the optimal m value corresponds to the peak of the correlation function, T, where PSFd1k and the scaled PSFd2k are the most relevant. However, based on the above formula, multiple variables are coupled to each other in near-field PSFs and we cannot derive the analytical solution of m with the peak of the correlation function between PSFd1k and PSFd2k. In this case, we simulated the PSFs by computer via Eq. (6) and search the most relevant areas by block matching to explore the statistical optimal solution of m.

The simulated PSFs are conducted as shown in Fig. 2, where the scattering process is created from a random phase mask based on the projection model [31,32]. A point light source was set at the on-axis position in the object plane as the corresponding point where Ok is located in the object plane. The simulated PSFs (normalized by removing the envelope) at different imaging distances through the modeled scattered layer are shown in Figs. 2(a) and 2(b). Visually, the intensity distribution of the two PSFs is strongly correlated, such as the green rectangles. Inspired by the wavefront slopes used from adaptive optics [33,34], we introduce the scaling vector into the algorithm to explore the detailed relationship between two imaging-distance PSFs, which is estimated by discrete block matching. We traverse PSFd2k to search the optimal area most relevant to a selected block in PSFd1k, and the translational relation from the selected block in PSFd1k to the optimal block in PSFd2k is equivalent to the scaling vector from the center point of the selected block in PSFd1k to PSFd2k. The block matching strategy can be expressed as where Nxd1,yd1 denotes the selected block with (xd1,yd1) as the center in PSFd1k; Mx,y denotes the searching block with (x,y) as the center in PSFd2k to match Nxd1,yd1; and the matched block in PSFd2k is centered around (xd2,yd2). Corr(·) calculates the correlation between two pixel blocks and we use the cross correlation in this paper [5]. In this way, a scaling vector is defined as an arrow with direction from (xd1,yd1) in PSFd1k to (xd2,yd2) in PSFd2k. To balance the computational complexity, we just utilized part of the discrete scaling vectors built by block matching, as shown in Fig. 2(c). The space between any two vectors vertically or horizontally is 20 pixels.

Then we calculate the m value of each estimated scaling vector in Fig. 2(c), which equals the distance between (xd2,yd2) and the scaling center divided by the distance between (xd1,yd1) and the scaling center, as shown in Eq. (7). Theoretically, all the scaling vectors which were estimated by block matching with maximum cross correlation share the same m value and the statistical analysis for Fig. 2(c) bears this out. The histogram distribution of m values as shown in Fig. 2(d) demonstrates that m stabilizes around 1.015 as a constant value, which meets the description of Eq. (7). Although the analytical solution of Eq. (7) is hard to derive, the scaling relationship between near-field PSFs statistical analysis does exist, which was proved by the above statistical analysis. We tentatively conclude that the on-axis (u=0,v=0) near-field scaling relationship between PSFd1k at Distance 1 plane and PSFd2k at Distance 2 plane can be expressed as where → denotes that PSFd1k is most relevant to the coordinate-scaled PSFd2k with the constant coordinate-scaled value m. Also, the on-axis scaling relationship in Eq. (9) can be generalized to the off-axis case as and the corresponding point where Ok is located changes from (u=0,v=0) to (u=uk,v=vk). As described in Eq. (10), for the off-axis scaling relationship situation, the scaling vector that starts at (x,y) in PSFd1k will end at [mx−(m−1)uk,my−(m−1)vk] in PSFd2k. This reveals that the scaling vectors have two features representing the location of Ok, which is also regarded as the scaling center from Distance 1 to Distance 2. First, the length of the scaling vector is proportional to the distance between the scaling vector and the scaling center. Second, the scaling center falls on the defined line by each scaling vector and each scaling vector points away from the scaling center (if d2>d1).

Under the spatially incoherent illumination, the speckle generated by one certain imaging target (supposed as Ok) is the linear superposition of many highly correlated PSFs. Actually, each PSF possesses a different scaling center corresponding to one point belonging to the imaging target. Considering the schematic of the proposed scattering imaging system in Fig. 1, the size of each imaging target falls within the ME range but the spacing between any two targets beyond. All the PSFs that form the single-target speckle possess the same scaling center approximately where the imaging target located around. Equation (11) explains the scaling relationship between the single-target near-field speckles as where (uk,vk) represents the scaling center of these two speckles which equals the location of Ok; (Δx,Δy) denotes the 2D coordinate differences between the scaling vector and the scaling center. Theoretically, the mentioned feathers between scaling vectors and scaling centers for PSFs still apply to the estimated scaling vectors by near-field speckles from (Δx+uk,Δy+vk) to (mΔx+uk,mΔy+vk).

After that, the scaling-vector-based detection algorithm, as shown in Fig. 3, is proposed based on Eq. (11) for multiple targets in Fig. 1 under the spatially incoherent light source. First, the scaling vectors are estimated by block matching as Eq. (8) from two imaging-distance multi-target speckles. We traverse Id2 to search the optimal area most relevant to a selected block in Id2. The density of estimated scaling vectors vertically or horizontally can be adjusted appropriately with the speckle resolution. Then, multiple targets mean that multiple scaling centers exist around the estimated scaling vectors. In the case that there is an uncertain number of scaling centers, we use the length information of each estimated scaling vector to determine some regions where the targets may locate. Any position whose scaling vector length is below a certain threshold is listed as the possible region. Next, the connected component analysis (8-connected) [35] will be applied for clustering these regions and the number of the connected components equals the number of imaging targets. The algorithm will choose the area with the minimum length of the scaling vector in each connected component as the rough location of each target. Finally, the direction information of the scaling vector would be used for each rough location to adjust the positioned row and column where the imaging target located accurately, because the line defined by each scaling vector belonging to one connected component theoretically passes through the target position in that component. The proposed scaling-vector-based detection algorithm is assisted with the block matching method and the connected component analysis to achieve blind target position detection, only exploiting two captured multi-target near-field speckles at different imaging distances. The detected position information includes the number, n, and the locations of the imaging targets corresponding to the object plane, (uk,vk),k=1,2,…,n, like the blue points in the Section 3.

After target position detection, the low-cross-talk region allocation strategy is proposed to extract the autocorrelation of each target, in order to simplify the infinite-solution problem in Eq. (3). We can select a small square region of the captured speckle, ε, in Id1 for autocorrelation calculation, which is centered by one detected target location (O1 as an example) with side length β as

Considering the envelope properties of speckle, the selected region affects the weight of each imaging target in Id1 by

This weight difference would be doubled if Id1 is transferred into the autocorrelation domain [25] as where Iε denotes the selected square region in Id1 belonging to ε; C(x−uk,y−vk) is spatially constant when ε is much smaller than the whole speckle to describe the remained weight for Ok caused by the envelope. Via the doubled weight difference, the autocorrelation of O1 can be extracted from the autocorrelation of the selected region in Id1 centered by the location of O1. Repeat the above steps n times for each target location and the autocorrelation of each imaging target will be reconstructed with low cross talk by other autocorrelations as shown in Eq. (14).

However, the extracted autocorrelations of each imaging target do carry some interference from other autocorrelations in theory, which leads to an unstable output after the traditional phase retrieval algorithm with the random phase as the initial input. In order to improve the stability of reconstruction, the modified phase retrieval algorithm is applied in the paper especially for the autocorrelation signal with interference. First, we blindly reconstruct a number of object images by the “hybrid input-output” and the “error-reduction” algorithms via different random initial phases [27]. Then, for each object image, the part whose intensity is less than 20% of the maximum intensity of that image is regarded as unstable noise and the intensity of that part is set to zero. Finally, the object image with the least change in the autocorrelation domain between the unprocessed one and the processed one is taken as the final optimal output.

Working with the modified phase retrieval algorithm, each target can be reconstructed with the help of the detected number and locations of the imaging targets. Then, these targets will be placed at the detected position to form a complete scene of the multiple targets exceeding the ME range without any prior information, eventually. It has been verified by the real experiments that the approximations that appear in this section are acceptable and the extracted position information is sufficient for a visually distinguishable reconstruction. In addition, we discuss the limitations of the proposed algorithm in the next section.

In the following section, we first describe the optical setup of the scattering imaging system, then detail the tests on a real scattering imaging system, and finally analyze the limitations of the proposed algorithm.

The multi-target large FOV scattering imaging system setup via the blind target position detection is shown in Fig. 4, which is extended from the single-shot scattering system proposed by Katz et al. [5]. A narrow bandwidth 532 nm single-frequency CW laser (Cobolt SambaTM-100) serves as the light source whose coherence is attenuated via rotating ground glass. One Thorlabs Optics 220-grit diffuser whose effective ME range is 16.6 mrad is placed between multiple targets and the sensor plane [25]. We used a single CMOS camera (Filr, pixel size=4.8μm, 1280×1024 pixels) on a 3D moving platform (DHC, minimum scale=10μm) to capture the large FOV speckles at different imaging distances.

First, various multiple targets were tested to demonstrate the effectiveness of the proposed algorithm. The left column shows the test on the mask “2FL” and the right column shows the test on a larger and more complex scene “01234.” Figures 5(a) and 5(f) describe the detailed parameters of multiple targets. The minimum distance between any two targets is 3.5 mm, and each target size is within 0.5 mm. In this part, the object distance (d0) equals 120 mm. The total scene size is much larger than the ME range (120mm×16.6mrad), but each target size is smaller than it. The real multi-target large FOV near-field speckles through scattering media at different imaging distances are captured via a 3D moving camera. Figures 5(c) and 5(h) are captured with d2=17.0mm and Figs. 5(d) and 5(i) with d1=16.5mm. Visually, no prior information (including the number and locations of the imaging targets) can be seen directly from these multi-target superimposed speckles. The red arrows in Fig. 5(e) show all the estimated scaling vectors using the block matching method based on the near-field speckles for mask “2FL,” which describe the scaling relationship between these two imaging-distance speckles. Obviously, there are three scaling centers that exist among these scaling vectors. Meanwhile, Fig. 5(e) shows the extracted connected components for mask “2FL” in the bottom right and the blue points denote the final detected multi-target locations. The number of the extracted connected components matches the number of the imaging targets and this multi-target scattering imaging problem is specified as the process of reconstructing three objects. After that, the autocorrelation of each imaging target as shown in Fig. 5(d) can be reconstructed via the operations as Eq. (14) in the selected speckle region centered by the detected position information in Fig. 5(e). Working with the modified phase retrieval algorithm, each imaging target can be reconstructed successfully and then put at its location to form a complete scene as shown in Fig. 5(b), which is visually close to the original one. The process of Figs. 5(h)–5(j) for target “01234” is the same as the process of Figs. 5(c)–5(e). The objective similarity evaluations between the reconstructions and the original targets by peak signal-to-noise ratio (PSNR) are shown in Table 1 and the detected target locations are shown in Figs. 5(b) and 5(g). Compared with the original scene in Figs. 5(a) and 5(f), these reconstructed results demonstrate that the multiple targets can be accurately positioned by scaling-vector-based position detection and the proposed method achieves multi-target large FOV blind reconstruction through scattering media in the real imaging system.

Second, to test the applicability of our proposed method for large FOV fluorescent biological observation through scattering media, the neuron-shape scattering imaging experiments are conducted. The neuron-shape mask scaled from the real dendrites of hippocampal neurons [36] has the complex shapes of targets and distribution irregularity of target locations, which increases the difficulty for reconstruction but satisfies the requirements for practical scattering scenes. The mask was set as the imaging targets in the proposed scattering imaging system with other experimental conditions consistent with Fig. 5. Figure 6(a) shows the neuron-shape mask for reconstruction and the detailed parameters. The reconstructed whole scene via the proposed method is shown in Fig. 6(b) and the main features are faithfully recovered. Certainly, no prior information can be seen directly from the captured large FOV speckle, as shown in Fig. 6(c), which was recorded with d1=16.5mm as the Distance 1 plane. In principle, the presented millimeter-scale experiments can be scaled to micrometer- or meter-scale scenarios for scattering imaging exceeding the ME range.

In the case that multiple targets are separated from each other as shown in Fig. 1, the proposed method blindly achieves multi-target localization and reconstruction exceeding the ME range. However, how to define the effective spacing between any two targets, beyond which the scaling-vector-based detection and scattering imaging algorithm works for multi-target speckles? To evaluate this limitation, we adjust the spacing between any two targets for mask “2FL” in Fig. 7(a) from 3.25 mm to 1.5 mm. Besides spacing, other parameters remain the same as those in Fig. 5(a), with the effective ME range keeping 2 mm in the object plane and the size of each imaging target falling within 0.5 mm. In this experimental environment, it is obvious that the reconstructed targets and locations are visually distinguishable when the spacing changes from 3.25 mm to 2.25 mm, as shown in Fig. 7(b). The estimated scaling vectors and detected locations with a spacing of 2.75 mm are shown in Fig. 7(d) as an example of reconstructions in good quality. As the spacing between any two targets reaches lower than 2.25 mm, the reconstructions are visually distorted, and the quality degrades. Meanwhile, the proposed blind position detection algorithm estimates the wrong number and the wrong locations of the imaging targets from the near-field speckles as shown in Fig. 7(e) with a spacing of 1.75 mm, which serves as misleading for the following reconstruction process. Objectively, the three-target averaged PSNRs curve between the reconstructions and the original targets with respect to the decreasing spacing is provided in Fig. 7(c). Theoretically, the autocorrelation of one selected target like Eq. (14) is noised more by other target autocorrelations with the decreasing spacing between any two targets, which makes the final reconstruction after phase retrieval worse. Meanwhile, according to Eq. (13), the captured multi-target speckle region centered by one detected target location (O1 as an example) is mainly composed of the speckle generated by O1 and the speckles generated by other imaging targets only make a small contribution for this region. Therefore, the scaling vectors in this region regard the location of O1 as the scaling center and there are multiple different scaling centers that exist between two imaging-distance speckles which correspond to the locations of different imaging targets. This is the reason why the proposed blind position detection algorithm can extract multi-target positions from two large FOV imaging-distance speckles, and Figs. 5(e) and 5(j) show the obvious boundaries of the speckle regions serving for different imaging targets. However, with the decreasing spacing, the weight differences in Eq. (13) are reduced and the scaling vectors in the speckle region centered by the location of O1 are noised more by other imaging targets. That will result in the low-accuracy target localization, and even the wrong number of the identified imaging targets.

As for the times that the proposed algorithm can exceed the ME range, it is not restricted by theory, only if the camera array can capture the full speckles when the number of imaging targets increases.

To summarize, we developed a multi-target large FOV scattering imaging method based on the blind target position detection. This technique only exploits two captured multi-target near-field speckles at different imaging distances, from which the target position cannot be seen directly. A major advantage of our approach is that target position information including the number and the locations of the imaging targets can be blindly reconstructed via the scaling-vector-based detection algorithm. After that, the autocorrelations can be calculated for speckle regions centered by the derived positions via the low-cross-talk region allocation strategy. Working with the modified phase retrieval algorithm, the whole scene of the multiple isolated targets exceeding the ME range can be successfully recovered. Unlike other methods of exceeding the ME range, no prior information of the target is required, which makes our technique more applicable for prior-information-free large FOV imaging and other scattering imaging techniques can cooperate with our proposed method to further improve the performance [25,26]. The real scattering imaging experiments demonstrate the effectiveness of the proposed method.

Actually, the detected target position information contains some errors compared with the initial multi-target positions. These errors come from three main sources: (1) the discrete sampling during camera acquisition and image processing; (2) the approximations of the theoretical derivation from Eq. (9) to Eq. (11) aiming to simplify the scattering model; and (3) a certain divergence angle of the experimental illumination, which makes the target positions change slightly from the object plane to the sensor plane. Actually, these errors did not affect the effectiveness of the proposed algorithm for target localization and reconstruction. On the other hand, the minimal effective spacing between any two targets for the proposed algorithm is not only limited by the target size and the properties of the scattering medium, but also the distance between the scattering medium and the camera sensor, and even the shape of the imaging target. In the future, a deeper research will be focused on this problem.

Additionally, the main time consumption of our proposed algorithm is spent in the process of estimating scaling vectors from two imaging-distance speckles and the time consumption will increase along with the increasing FOV of the captured speckles. When the resolution of the captured speckles in Figs. 5(c) and 5(d) is 1600×2500 pixels, the length of the selected square region (β) is set as 400 pixels and the total time consumption in MATLAB 2018b is about 6.1 h. In addition, in the image processing part, the raw captured multi-target speckle with the slowly varying envelope is normalized by dividing the raw speckle by a low-pass-filtered version of it before calculating autocorrelations [5]. Especially in the process of estimating scaling vectors, the normalized speckle is further smoothed by a low-pass filter to remove the existing noises by camera which are spatially invariant at different imaging distances and improve the accuracy for localization.

Finally, the imaging distances, d1 and d2, are not fixed in the proposed method, so as the gap between two near-field speckles, d2−d1. In addition, these distance parameters were not used as the known information for reconstruction in this paper. In practice, the imaging distance can be adjusted to satisfy the requirements for different scenes, as long as the near-field scaling relationship between two imaging-distance speckles exists, and the light field techniques [37] may be applied to the proposed system to speed up the acquisition process. Meanwhile, the relationship between parameter m and multiple variables (including d0, d1, d2, and TM) via statistical analysis will be the research plan in the future, aiming to dig out more useful information from two near-field speckles at different imaging distances. Furthermore, the proposed method can still work theoretically when the imaging targets are sandwiched between two scattering layers or looked around corners [5], which will be tested experimentally in the future work, and the raw data used to generate the results presented in this manuscript is available at https://cloud.tsinghua.edu.cn/d/296c066dbcc243839f52/.

Acknowledgment. We would like to thank Xiangsheng Xie for helpful discussions.