Non-invasive imaging through strongly scattering media based on speckle pattern estimation and deconvolution

Imaging through scattering media is still a formidable challenge with widespread applications ranging from biomedical imaging to remote sensing. Recent research progresses provide several feasible solutions, which are hampered by limited complexity of targets, invasiveness of data collection process and lack of robustness for reconstruction. In this paper, we show that the complex to-be-observed targets can be non-invasively reconstructed with fine details. Training targets, which can be directly reconstructed by speckle correlation and phase retrieval, are utilized as the input of the proposed speckle pattern estimation model, in which speckle modeling and constrained least square optimization are applied to estimate the distribution of the speckle pattern. Reconstructions for to-be-observed targets are realized by deconvoluting the estimated speckle pattern from the acquired integrated intensity matrices (IIMs). The qualities of reconstructed results are ensured by the stable statistical property and memory effect of laser speckle patterns. Experimental results show that the proposed method can reconstruct complex targets in high quality and the reconstruction performance is robust even much less data are acquired.

while the main problem is how to acquire the point-spread speckle pattern. Point-spread-function (PSF) capturing based deconvolution method 19,20 needs to put a point source behind the scattering media during the estimation process, whose invasiveness impedes the practicability for many real applications.
To non-invasively image complex targets through strongly scattering media with resolution beyond the limitation of illumination speckle grain lighted on the targets, a preliminary speckle pattern estimation and deconvolution method was proposed in our previous work 47 . The method is extended in this paper by adding cross-correlation matching method for restoring orientations and positions, deriving theoretical model for the speckle pattern estimation, analyzing the solution space and data reduction and comparing experimental results for imaging targets with different sizes and structural complexities. Imaging system based on speckle correlation 15 is utilized to collect IIMs of the to-be-observed targets. Exploiting the frequency-bounded characteristics of the speckle pattern 48 and the convolutional relationship between the imaging target and the speckle pattern, speckle estimation based on the constrained least square optimization and speckle modeling is proposed to estimate the speckle pattern distribution. IIMs are aligned by the orientational and positional information recovered from the cross-correlation matching method. The complexity in solving the speckle pattern estimation model and the quality of estimation using the training targets in different shapes with different combinations are discussed in detail. Deconvolution is performed to reconstruct targets using their corresponding IIMs and the estimated speckle pattern. Reconstruction using much less acquired IIM data is also analyzed to demonstrate the feasibility in reducing the data collection complexity. Compared with speckle correlation [15][16][17] and bispectrum analysis 18 based approaches, in which the stochastic perturbations in point-spread speckle pattern are considered as statistic noises, our method directly estimates the speckle pattern and reconstruct the target by deconvolution algorithm. The demand of high-resolution angular sampling is relieved and it can recover high-quality details of complex targets with directional information. Experimental results and discussions on real-captured data are provided to verify the effectiveness of the proposed method. Compared with the existing approaches, the proposed method provides much higher reconstruction quality to much more complex to-be-observed targets with much larger FoV.

Principles and Methods
Experimental setup and denotations. The schematic diagram of data collection process is illustrated in Fig. 1. In the enclosed environment surrounded by the scattering media, several imaging targets with different structural complexities may appear in the imaging region. The targets that can be reconstructed by phase-retrieval algorithm, namely training targets, are denoted by . And the to-be-observed targets are denoted by O c . Imaging system based on raster scanning and speckle correlation, which was first proposed by Bertolotti et al. 15 , is utilized to collect the IIMs for the imaging targets. When imaging target is in the common imaging region behind the scattering media, the laser beam with incident angle of θ θ ( , ) x y will generate speckle pattern lighting on the target through the scattering media. In case of the existence of memory effect [21][22][23] , the point-spread speckle pattern can be approximately expressed by , ) x y 1 1 , where d 1 is the distance between the imaging target and the scattering media. The intensity of the speckle pattern transmitted and diffused by the scattering media, which encodes the information of the imaging target, is recorded by the photon detector. Imaging target can be generally denoted by ( 2) , where (x, y) is the coordinate on the object plane and r r ( , ) (1) (2) is its spatial position relative to the center of the FoV. The recorded integrated intensity corresponding to a specific incident angle of laser beam, θ θ ( , ) x y , can be expressed as x y x y x y (1) ( 2) 1 1 1 (1) 1 (2) where * denotes the convolution operator and IIM is the integrated intensity matrix. The pixel intensity of IIM corresponds to the total intensity of the image captured at a laser incident angle. The sample IIMs collected for O s Fig. 1(a-c) are shown in Fig. 1(d-f), respectively.
Speckle pattern estimation. From Eq. (1), it can be observed that once the point-spread speckle pattern that scans the imaging region is known, theoretically, the target can be reconstructed by deconvoluting S from IIM by 47 where e i represents the error caused by imaging noise and phase-retrieval algorithm.
x y is the dimension of the imaging target and × K K x y is the number of sampling points in IIM s i . It is a linear system with  (3) is an underdetermined problem and e i is also unknown. Considering the intensity distribution of speckle pattern is always positive and it presents frequency-bounded characteristics inherited from Fourier transformation of the laser beam wavefront 48 , the prior of speckle pattern distribution is introduced to build a constrained least square model for speckle estimation as 47 where ∆ constrains the variation between the intensity of two adjacent points in the speckle pattern. ∆ is inversely proportional to the size of the illumination speckle grain, and the speckle grain size is proportional to d 1 and wavelength of the laser while inversely proportional to the diameter of laser beam 48 . The speckle estimation model in Eq. (4) is a convex problem 49 , which can be solved by the solvers of convex optimization problem (CVX 50,51 is used in our experiments).

Phase retrieval and cross-correlation matching for reconstruction of training targets. The
phase-retrieval reconstruction of training target, denoted by O s i , is used as the input of Eq. (4). To reconstruct O s i , the sharply peaked property of the average autocorrelation of speckle pattern 52 is exploited first. Through autocorrelation of IIM and plugging Eq. (1) into it, the Fourier magnitude spectrum of O s i is given by 15  where * denotes the correlation operator. Then, phase-retrieval algorithm, Hybrid Input-Output and Error-Reduction 15,45,46 , is applied to recover O s i via iteratively updating the phases with fixed magnitude spectrum.
After phase retrieval, the orientational and positional information of O s i is lost. Also, the training targets may not be spatially aligned in the imaging region. It leads to the spatial disconformity between the captured IIMs with the convolutions of O s i and S. Directly putting O s i and the captured IIM s i into Eq. (4) will result in mismatches and cannot estimate the local spatial distributions of speckle pattern accurately. Thus, cross-correlation matching is proposed to correct the orientations of O s i and to align IIM s i in Eq. (4). According to Eq. (1), the IIM of O s i can be expressed as is the spatial position of O s i relative to the center of FoV. The autocorrelation of speckle pattern is a sharply peaked function, which can be expressed by 48 where δ denotes impulse function; B and C are the constants depending on the size of speckle grain. Applying cross-correlation to the IIMs of two training targets O s i and O s j , we have Thus, by comparing the four possible orientation combinations of O s i and O s j , the orientation which maximizes the correlation between the correlation of rotated phase-retrieved results and the correlation of IIMs is achieved by As a result, the reconstructed training targets with correct orientations are obtained by i is an approximation of O s i with some reconstruction errors, the spatial relative offsets between O s i and O s j can be estimated by Finally, the speckle estimation model in Eq. (4)

Model Analysis
Training targets analysis. The speckle estimation performance of Eq. (15) is decided by the shape, size and the number of training targets. The increase in the number of training targets, the larger shape differences among training targets and the size reduction of each target will result in better estimation performance. However, the estimation qualities of entries in matrix S are not the same. The objective function of the estimation model, Eq. (15), is the least square solution of linear equations, s i s i in which the entries of S in the central are included in more independent equations than those in the margins. Thus, the central entries have smaller solution spaces and higher probabilities to be solved with lower errors. Eq. (16) can also be discretely formed by where χ θ θ y y x y x x y y x y y y y s i x y x y x y 1 1 x y , as shown in Fig. 2(a-d), are selected. Speckle patterns at object plane are generated by simulating the wave propagation of random phased laser beam with Fresnel diffraction equation 43 . Raster scanning process (corresponds to changing the incident angle of laser beam with the memory effect 21-23 ) is simulated to acquire the IIMs in size of 60 × 60 = K K ( , 60) x y . To get rid of the reconstruction error caused by phase-retrieval algorithm, original training targets are utilized to estimate the speckle patterns using the proposed model in Eq. (15). Cropping the estimated speckle patterns in size of 60 × 60 in the center, i.e. the region lined in white, the peak signal to noise ratio (PSNR) is calculated using the simulated ground truth, as an example shown in Fig. 2(e), as the reference. 100 simulated speckle patterns are tested to average the PSNR for a robust result. As we can see from the results, the estimation quality is positively related to the number of the training targets used. For example, the estimation quality using two training targets is much higher than only using one training target, while much lower than using three or four training targets. Also, it can be found that the estimation quality is mainly determined by the rank-variable ratio, i.e. χ = rank A NoV A b ( )/ ( ). When the number of training targets is increased, the increase in the number of independent equations, i.e. rank(A), is larger than that in the number of variables, i.e.
) . And it results in the increase of rank-variable ratio. Since the higher rank-variable ratio corresponds to the less underdetermined variables, the estimation quality can be improved obviously. The variation in the shape of the training targets also results in the change of the rank-variable ratio. If the targets with higher structural complexities can be reconstructed well by phase retrieval, the bigger structural differences among the training targets result in the higher rank(A), which benefits the speckle pattern estimation.
Reconstruction with less data. After getting the estimated speckle pattern S from Eq. (15), deconvolution can be executed for the acquired IIM to reconstruct an imaging target with much higher accuracy and robustness than the phase-retrieval methods. It is also observed that we can further reduce the data acquisition complexity θ θ  ( 2) and S x y ( , ) with the sampling intervals as θ θ ∆ ∆ ( , ) , ) Thus, using the down-sampled estimated speckle pattern, denoted by S , it is still possible to reconstruct the imaging target in lower resolution by deconvoluting IIM c with S .

Results and Discussion
In this section, the performance of the proposed speckle pattern estimation model is demonstrated by testing on a real imaging system, in which two Edmund Optics 120-grit ground-glass diffusers 15 are used as the scattering media. In this system, speckle scanning and recovery method proposed by Bertolotti et al. 15 is used. The laser beam is controlled by two galvanometric scanners in vertical and horizontal directions separately to generate different incident angles at the scattering media. The details of the imaging system can be found in the electronic supplementary material. To reduce the noises and artifacts, Richardson-Lucy deconvolution algorithm 53,54 is applied for reconstruction.
Four training targets, as Fig. 3(a), are used in the experiments. Their sizes range from 20 microns to 25 microns. The distances between the targets and scattering media, d 1 in Eq. (1), are set as 5 millimeters. Scanning angles are divided evenly into 600 intervals, ranging from − .°2 1 to .°2 1 in both horizontal and vertical directions. Thus, the corresponding IIMs collected are shown in Fig. 3 to be those in Fig. 3(e) to reduce the computational complexity of solving Eq. (15). Figure 3(f) and (g) show the estimated speckle patterns using only two training targets and four training targets respectively. The ground truth of the speckle pattern is shown in Fig. 3(h). Comparing Fig. 3(f) and 3(g) with the ground truth, it can be found that both the two estimated speckle patterns are visually similar to the ground truth especially considering the low-frequency features. The result estimated using four training targets preserves more high-frequency components, which may lead to reconstructing more details of the imaging targets by deconvolution.
Applying the estimated speckle pattern to the deconvolution process, the to-be-observed targets with complex structures can be reconstructed. The reconstruction quality is compared with those generated by phase retrieval 15,45,46 to demonstrate the efficiency of the proposed method. As shown in Fig. 4(a), eight complex targets with different textural complexities and different number of subsections are tested. Their sizes range from 40 microns to 165 microns. IIMs of these targets are collected in the same way with the training targets. The reconstructed results are compared visually with the original images. Also, structural similarity index (SSIM) 55 is calculated between the reconstructed results and the original images to evaluate the structural reconstruction capability objectively. SSIM closer to 1 corresponds to a higher reconstruction quality. Figure 4(b) shows the reconstructed results using phase retrieval, each of which is the best result selected from performing 20 times of phase retrieval using random initial phases. It can be found that the quality of the results deteriorates sharply when targets become complex or have separated subsections. Also, because of the randomness of initial phases, the positions and shapes of the reconstructed results change a lot for each test. Figure 4(c,d) show the reconstructed results generated by deconvoluting estimated speckle pattern in Fig. 3(f) and that in Fig. 3(g), respectively. It is obvious to see that compared with the results generated by phase retrieval, the proposed method improves the reconstruction performance for all the targets significantly, which is also robust to the structure of the targets. Although the reconstructed results derived by deconvoluting speckle pattern estimated from two training targets lose some high-frequency details, as shown in Fig. 4(c), their reconstruction qualities are still much higher than that of phase retrieval. It demonstrates that an acceptable reconstruction quality can be achieved by applying less training targets. For the targets with separated subsections, like those on the last two rows with spatial extensions more than 150 microns, the proposed method always provides reconstruction results superior to those of phase retrieval. It demonstrates a big improvement in the FoV for our method. Figure 4(e-g) show the reconstruction results generated by deconvolution using 4, 16, 36 times lower sampling rate in collecting the IIMs of the imaging targets, as defined in Eq. (22). The speckle pattern estimated using four training targets, as shown in Fig. 3(g), is also down-sampled to be S . The reconstructed results are up-sampled using bilinear interpolation and SSIM is calculated between the up-sampled results and the original images. It can be found that acquiring only 1/4 × 1/4 of the data, the reconstruction quality of the proposed method is still acceptable. Even if reducing the data acquired to be only 1/6 × 1/6 of those used by phase retrieval, the proposed Scientific RePoRtS | (2018) 8:9088 | DOI:10.1038/s41598-018-27467-1 method can still provide reconstruction quality better than that of phase retrieval. It shows high potential in significantly reducing the complexity in speckle scanning and data acquisition process for collecting IIMs, which greatly benefits the real applications.

Conclusions
To non-invasively image complex targets through strongly scattering media with resolution beyond the limitation of illumination speckle grain lighted on the targets, speckle pattern estimation and deconvolution method is proposed. Cross-correlation matching is proposed to correct the orientations of training targets and align their corresponding IIMs, which are applied as the inputs of the speckle pattern estimation model. Both theoretical derivation and experimental results have demonstrated that the proposed speckle pattern estimation model can accurately estimate the distribution of the point-spread speckle pattern. The proposed deconvolution method is also robust even with much less data acquired. In conclusion, our imaging method relieves the complexity limitation for imaging targets, improves the reconstruction performance and reduces the data-collection complexity for the to-be-observed targets.
The principal limitation of our proposed work is the dependency on the reconstruction accuracy of training targets by phase retrieval. While, the reconstruction accuracy of phase retrieval is constrained by the relative spectral bandwidths between the imaging targets and the speckle pattern 17 . So, we put this as one of our future works to seek the way, such as introducing low rank feature of S, to improve the prior information for the speckle pattern estimation model, and to further analyze the solution space. Evaluating structural complexity of target from its IIM can be quite valuable in real scenario for selecting suitable training targets without phase retrieval. A preliminary experiment shows a kind of correlation between the entropy of the IIM and the phase-retrieval performance, which will be investigated further with other possible features. Another problem faced in the data collection process is the movements of imaging targets when collecting their IIMs, and it causes the structural deformation of IIMs and reconstructed results. A possible solution can be used is to track targets in real-time by analyzing the linear relationships between the variance and decorrelation time of the integrated intensity with target displacement along the axial and transversal directions, respectively 56 . Thus, the incident-angle offset of laser beam can be calculated. Beyond these, deconvoluting to-be-observed targets from different depths and introducing the deformation of memory effect can be future works as well.  Fig. 3(g).